#整體設定，含載入套件
source("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/R4BS_setup.R")

#讀取檔案
dta <- read.csv('../Data/KIT_memory.csv', 
                header=T, stringsAsFactors = TRUE)

#看資料結構
#程式報表10.1
dim(dta)

#檢視資料格式與前六筆
head(dta)

#寬資料轉成長資料
#也另外製造「月」，以及將月置中的「月置中」變項
dtaL <- dta |> 
  tidyr::pivot_longer(cols=3:7,
                      names_to ='月齡',  
                      values_to = '記憶分數') |>
  dplyr::mutate(月 = parse_number(月齡), 
                月置中 = scale(月, scale = F)[,1])

#程式報表10.2
head(dtaL)

#看五波的平均數與變異數，母親教育程度的影響在月齡較大時顯現
#程式報表10.3
dtaL |> 
  tidyr::unite("月齡母親教育程度", 月齡, 母親教育程度) |>
  dplyr::select("月齡母親教育程度", "記憶分數") |>
  gtsummary::tbl_summary(by=月齡母親教育程度,
            digits=list(記憶分數 ~ c(1, 2)),
            type = all_continuous() ~ "continuous2",
            statistic = list(all_continuous() ~ "{mean} ({sd})")) |>                          
  modify_header(label ~ "**變項**") 

#看看偏態與峰度
#程式報表10.4前
dtaL |> 
  dplyr::group_by(月齡, 母親教育程度) |>
  dplyr::reframe(平均 = mean(記憶分數),
                 標準差 = sd(記憶分數),
                 偏度 = moments::skewness(記憶分數),
                 峰度 = moments::kurtosis(記憶分數))

#歷次分數間相關
#程式報表10.4後
dta |> 
  dplyr::select(where(is.numeric)) |> 
  cor() |> 
  round(3)

#繪製不同母親教育程度五波圖
#圖10.1a
p1 <- ggplot(data = dtaL, 
       aes(x = 月, y = 記憶分數)) +
 ggbeeswarm::geom_quasirandom(aes(group=月), alpha=.1)+
 facet_wrap(vars(母親教育程度))+
 stat_summary(fun = mean, geom = 'line', linewidth=.8) +
 scale_x_continuous(breaks=c(3,6,12,18,24))+
 labs(x = '月齡', 
      y = '記憶能力分數') +
 theme(legend.position = 'top')

#看五波的機率密度函數圖
#圖10.1b
p2 <- ggplot(dtaL, 
       aes(x=記憶分數, y=月齡)) +
  ggridges::geom_density_ridges2()+
  scale_x_continuous(breaks=seq(10, 70, by=10))+
  facet_wrap(vars(母親教育程度))+
  labs(x = '記憶能力分數', 
       y = '月齡')
p1+p2

#不同母親教育程度下的實徵累積機率密度圖
#圖10.2
ggplot(dtaL,
       aes(x=記憶分數, group=母親教育程度)) +
  stat_ecdf()+
  facet_wrap(vars(月齡))+
  labs(x="記憶能力分數",
       y="實徵累積機率密度")

#看看個別資料，畫上迴歸線並配上區間
#圖10.3
ggplot(data = dtaL, 
       aes(y=記憶分數, x=月)) +
 geom_line(aes(group = 識別碼), 
           linewidth = .2,
           alpha=.5,
           col=8) +
 stat_summary(fun.data = 'mean_cl_boot', geom = "errorbar", width = 0.5) +
 stat_smooth(method="lm", 
             formula= y ~ poly(x, 2),
             se=FALSE, 
             col=1, 
             linewidth=.5)+
 facet_wrap(vars(母親教育程度))+
 scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
 scale_y_continuous(breaks=seq(8, 72, by=8))+
 labs(x = '月齡', 
      y = '記憶能力分數')

#前述圖形看起來不完全是線性，後續進一步探索
#將資料框轉成 tsibble，方便後續使用
dta_tb <- brolgar::as_tsibble(x = dtaL,
                             key = 識別碼,
                             index = 月,
                             regular = FALSE)

#brolgar套件的features功能，可以計算每個幼兒記憶分數的五數綜合（含中位數）
#利用gghighlight，在繪圖時，特別強調中位數特別高或低的幼兒資料
#圖10.4
dta_tb |>
  brolgar::features(記憶分數, feat_five_num) |>
  dplyr::left_join(dta_tb, by = '識別碼') |>
  ggplot(aes(x = 月,
             y = 記憶分數,
             group = 識別碼)) +
  geom_line(linewidth=.2) + 
  gghighlight::gghighlight(med > 50 | med < 28, 
                           label_key = 識別碼,
                           use_direct_label = FALSE)+
   scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
  labs(x = '月齡', 
      y = '記憶能力分數')

#brolgar套件中的facet_strata()，預設將資料分成12組，助於探索資料
#看看個別資料
#圖10.5
dta_tb |>
  ggplot(aes(x = 月,
             y = 記憶分數,
             group = 識別碼)) + 
  #geom_point()+
  geom_line(alpha=.2, linewidth=.2)+
  brolgar::facet_strata()+
  scale_x_continuous(breaks=c(3, 6, 12, 18, 24))

#隨機取23位來看一下資料
#前面觀察到資料並非線性，嘗試以二次式迴歸配適資料
#為對比二次式，也用樣條迴歸（spline regression）配適資料
#圖10.6
set.seed(18072023)
dta |> dplyr::slice_sample(n=23) |>
  tidyr::pivot_longer(cols=3:7,
                      names_to ='月齡',  
                      values_to = '記憶分數') |>
  dplyr::mutate(月 = parse_number(月齡)) |>
  ggplot()+
  aes(x=月, y=記憶分數)+
  geom_point(pch=1, size=rel(3))+
  stat_smooth(method='lm', formula= y ~ poly(x, 2), 
              col=8,
              alpha=.5,
              linewidth=.6,
              se=FALSE)+
  stat_smooth(method='lm', 
         formula = y ~ splines::bs(x, knots = c(3, 12, 24), degree = 1),
              col=2,
              linetype='dotted',
              linewidth=1,
              se=FALSE)+
  facet_wrap(vars(識別碼)) +
  scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
  scale_y_continuous(breaks=seq(16, 72, by=16))+
  labs(x = '月齡', 
       y = '記憶能力分數',
       title="實灰:二次多項式, 點紅:樣條線性迴歸")

#以nlme，做每個人的迴歸線
m1 <- nlme::lmList(記憶分數 ~ 月置中 + I(月置中^2) | 識別碼, data = dtaL, subset=母親教育程度=='專科以下') |> 
  coef() |> 
  as.data.frame() |>
  dplyr::mutate(母親教育程度='專科以下')

m2 <- nlme::lmList(記憶分數 ~ 月置中 + I(月置中^2) | 識別碼, data = dtaL, subset=母親教育程度=='大學以上') |> 
  coef() |> 
  as.data.frame() |>
  dplyr::mutate(母親教育程度='大學以上')

m12 <- rbind(m1, m2)
names(m12)[1:3] <- c("b0", "b1", "b2")

#截距、一次項與二次項係數間的相關
cor(m12[,-4])

#看看母親教育程度對係數的影響
#程式報表10.5
m12 |> 
  gtsummary::tbl_summary(by=母親教育程度,
                         statistic = list(all_continuous() ~ "{mean} ({sd})"))

#不同母親教育程度所得之截距，線性斜率，二次項效果估計值
#圖10.7
g1<-ggplot(m12,
       aes(x=b0, y=b1, color=母親教育程度))+
  geom_point(pch=1)+
  stat_ellipse()+
  scale_color_grey()+
  labs(x='截距估計值', 
       y = '線性斜率估計值')+
  theme(legend.position='top')
g2<-ggplot(m12,
       aes(x=b1, y=b2, color=母親教育程度))+
  geom_point(pch=1)+
  stat_ellipse()+
  scale_color_grey()+
  labs(x='線性斜率估計值', 
       y = '二次項效果估計值')+
  theme(legend.position='top')
g3 <- ggplot(m12,
       aes(x=b0, y=b2, color=母親教育程度))+
  geom_point(pch=1)+
  stat_ellipse()+
  scale_color_grey()+
  labs(x='截距估計值', 
       y = '二次項效果估計值')+
  theme(legend.position='top')
g1 + g2 / g3

#潛在成長模型可以視為因素分析的子模型，我們以 lavaan 進行分析
#先選擇母親教育程度為大學以上者分析
dta_2 <- dta |> dplyr::filter(母親教育程度 == '大學以上')

#先試試看線性模型
growth1 <- '
intercept =~ 1*月_03 + 1*月_06 + 1*月_12 + 1*月_18 + 1*月_24
slope =~ 0*月_03 + 1*月_06 + 3*月_12 + 5*月_18 + 7*月_24
'
#配適模型
fit1 <- lavaan::growth(model = growth1, data = dta_2)
#程式報表10.6前，後
fit1 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))

fit1 |> parameterestimates() 

#試試看二次模型

growth2 <- '
intercept =~ 1*月_03 + 1*月_06    + 1*月_12 + 1*月_18 +  1*月_24
linear =~ (-3)*月_03 + (-2)*月_06 + 0*月_12 + 2*月_18 +  4*月_24
qudratic =~  9*月_03 + 4*月_06    + 0*月_12 + 4*月_18 + 16*月_24
'
fit2 <- growth(model = growth2, data = dta_2)

#程式報表10.7前
fit2 |> parameterestimates() |>
  filter(!is.na(z)) 

#修改二次模型，要求二次參數變異數為小的正數
#程式報表10.10後
growth21 <- '
intercept =~ 1*月_03 + 1*月_06    + 1*月_12 + 1*月_18 +  1*月_24
linear =~ (-3)*月_03 + (-2)*月_06 + 0*月_12 + 2*月_18 +  4*月_24
qudratic =~  9*月_03 + 4*月_06    + 0*月_12 + 4*月_18 + 16*月_24
qudratic~~.01*qudratic
'
fit21 <- growth(model = growth21, data = dta_2)
#程式報表10.7中，後
fit21 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))

fit21 |> parameterestimates() |>
  filter(!is.na(z)) 

#線性模型與二次模型比較
#程式報表10.8
anova(fit1, fit21)

#包含形狀因素的潛在成長模型，形狀不固定
growth3 <- '
intercept =~ 1*月_03 + 1*月_06 + 1*月_12 + 1*月_18 + 1*月_24
shape =~     0*月_03 + 月_06   +   月_12 + 月_18   + 1*月_24
'
fit3 <- growth(model = growth3, data = dta_2)
#程式報表10.9前，後
fit3 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))

fit3 |> parameterestimates() |>
  filter(!is.na(z)) 

#程式報表10.10
semTools::compareFit(linear = fit1, quadratic = fit21, shape = fit3, nested=F) |> summary()

#繪製路徑圖
#將圖存為 .png
#圖10.8
pl_01 <- lavaanPlot(model=fit3,
           edge_options = list(color = "grey"),
           coefs = TRUE,
           stand = FALSE,
           covs = TRUE)
dir.create(file.path(here::here()), 'figure')
lavaanPlot::save_png(pl_01, "figure/plot_01.png")

knitr::include_graphics("figure/plot_01.png")

#我們選擇包含形狀因素的潛在成長模型
#進一步看看母親教育程度對截距與形狀的影響。
#先看看不同母親教育程度下，資料配適情形。
dta_1 <- dta |> 
  dplyr::filter(母親教育程度 == '專科以下')

fit32 <- growth(model = growth3, data=dta_1)
#程式報表10.11前，後
fit32 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))

fit32 |> parameterestimates() |>
  filter(!is.na(z)) 

#將母親教育程度變成 0,1 的虛擬變項
dta <- dta |> 
  dplyr::mutate(母親大學以上 = ifelse(母親教育程度 =='大學以上', 1 , 0))

#包含形狀因素的潛在成長模型，讓性別影響截距與形狀因素。
growth4 <- '
intercept =~ 1*月_03 + 1*月_06 + 1*月_12 + 1*月_18 + 1*月_24
shape =~     0*月_03 + 月_06   + 月_12  + 月_18   + 1*月_24
intercept ~ 母親大學以上
shape ~ 母親大學以上
'
fit4 <- growth(model = growth4, data = dta, fixed.x = T)
#程式報表10.12前，後
fit4 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))

fit4 |> parameterestimates() |>
  filter(!is.na(z)) 

#比較實際資料與預測值平均
#計算實際平均數
df <- aggregate(cbind(月_03, 月_06, 月_12, 月_18, 月_24) ~ 母親教育程度, 
                data = dta, mean) |>
  tidyr::pivot_longer(cols=-1, 
                      names_to = "月", 
                      values_to = '記憶分數平均') |>
  dplyr::mutate(月 = parse_number(月))
#計算預測值
t <- c(0,.145, .618, .808, 1)
df$估計值平均 <- c((20.506-.133)+(32.870+2.262)*t, 20.506+32.870*t)

#畫圖
#圖10.9
df |> tidyr::pivot_longer(cols=3:4,
                      names_to = "估計",
                      values_to = "平均分數") |>
 ggplot()+
   aes(y=平均分數, x=月, group=估計, color=估計)+
   geom_point(size=rel(5), pch=1)+
   geom_line(linewidth=.4)+
   facet_wrap(vars(母親教育程度))+
   scale_color_grey(end=.6)+
   scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
   labs(x = '月齡',
        y = '記憶能力平均分數')+
   theme(legend.position='top')

#繪製路徑圖
#將圖存為 .png
#圖10.10
pl_02 <- lavaanPlot(model=fit4,
           edge_options = list(color = "grey"),
           coefs = TRUE,
           stand = FALSE,
           covs = TRUE)
dir.create(file.path(here::here()), 'figure')
lavaanPlot::save_png(pl_02, "figure/plot_02.png")

knitr::include_graphics("figure/plot_02.png")